Magnetization Switching in Nanowires: Monte Carlo Study with Fast Fourier 

Transformation for Dipolar Fields 
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For the investigations of thermally activated magnetization reversal in systems of classical mag- 
netic moments numerical methods are desirable. We present numerical studies which base on time 
quantified Monte Carlo methods where the long-range dipole-dipole interaction is calculated with 
the aid of fast Fourier transformation. As an example, we study models for ferromagnetic nanowires 
comparing our numerical results for the characteristic time of the reversal process also with numer- 
ical data from Langevin dynamics simulations where the fast Fourier transformation method is well 
established. Depending on the system geometry different reversal mechanism occur like coherent 
rotation, nucleation, and curling. 
Pacs: 75.10.Hk, 75.40.Mg, 75.40.Gb 
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I. INTRODUCTION 

Small magnetic particles in the nanometer regime are 
interesting for fundamental research as well as for appli- 
cations in magnetic devices. With decreasing size ther- 
mal activation becomes more and more important for the 
stability of nanoparticles and is hence investigated ex- 
perimentally as well as theoretically. Wernsdorfer et al. 
studied magnetization reversal in nanowires 0] as well as 
in nanoparticles j2| experimentally. They found that for 
very small particles the magnetic moments rotate coher- 
ently as in the Stoner-Wohlfarth model g while for larger 
system sizes more complicated nucleation processes oc- 
cur. 

Numerical studies of the thermal activation in mag- 
netic systems base either on Langevin dynamics J^-q] 
or on Monte Carlo methods. With the second method, 
mainly nucleation processes in Ising systems have been 
studied but also vector spin models have been used 
to investigate the switching behavior in systems with con- 
tinuous degrees of freedom |p|-|l3|. In these numerical 
studies the dipole-dipole interaction is often neglected or 
at least approximated due to its large computational ef- 
fort. 

In this work, we will use Monte Carlo methods in or- 
der to investigate the thermally activated magnetization 
behavior of a systems of magnetic moments including the 
dipole-dipole interactions. Since the calculation of these 
long-range interaction is extremely time consuming when 
performed straightforward it is necessary to develop ef- 
ficient numerical methods for this task. We will demon- 
strate that the implementation of fast Fourier transfor- 
mation (FFT) methods — which are already established 
in the context of micromagnetic simulations based on the 
Landau-Lifshitz equation of motion Jl4|,[l5| — is possible 
also in a Monte Carlo algorithm where a single-spin-flip 
method is used. 

Our approach is applied to the important problems 



of thermally activated magnetization reversal in a model 
for nanowires. Depending on the system geometry, mate- 
rial parameters, and the magnetic field different reversal 
mechanisms occur. Very small wires reverse by a coher- 
ent rotation mode while for sufficiently long wires the 
reversal is dominated by nucleation With in- 

creasing system width an additional crossover sets in to 
a reversal by a curling-like mode. Our numerical results 
based on Monte Carlo simulations are supplemented by 
Langevin dynamics simulations for comparison. 



II. MODEL FOR A NANOWIRE 

We consider a classical Heisenberg Hamiltonian for lo- 
calized spins on a lattice, 

H = -J^S^S J -/i s B-^S i 

(ij) i 



(1) 



where the Sj = /Ltj/^s are three dimensional magnetic 
moments of unit length. The first sum represents the 
ferromagnetic exchange of the moments where J is the 
coupling constant, the second sum is the coupling of 
the magnetic moments to an external magnetic field 
B, and the last sum represents the dipolar interaction. 
w = /ioMs/(4™ 3 ) describes the strength of the dipole- 
dipole interaction with no = 4tt ■ 10 _7 Vs(Am) -1 . We 
consider a cubic lattice with lattice constant a. The 
are unit vectors pointing from lattice site i to j. 

We model a finite, cylindric system with diameter D 
and length L (in number of lattice sites). All our simula- 
tions start with a configuration where all spins point up 
(into z direction) . The external magnetic field is antipar- 
allel to the magnetic moments, hence favoring a reversal 
process. As in earlier publications Jl0| |l3[| we measure the 
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characteristic time of the reversal process, i. e. the time 
when the z component of the magnetization of the sys- 
tem changes its sign averaged over either 100 (Langevin 
dynamics) or 1000 (Monte Carlo) simulation runs. 



III. NUMERICAL METHODS 

In the following we will mainly use a Monte Carlo 
method with a time step quantification as derived in |l2j ] 
and later further discussed in fll3| . The method is based 
on a single-spin-flip algorithm. The trial step of the al- 
gorithm is a random movement of the magnetic moment 
within a certain maximum angle. In order to achieve 
this efficiently we construct a random vector with con- 
stant probability distribution within a sphere of radius R 
and add this vector to the initial moment. Subsequently 
the resulting vector is normalized. The trial step width R 
is connected to a time scale At according to the relation 



r2 = JOksTar, At 



(l + a 2 K 



(2) 



Using this relation one Monte Carlo step per spin (MCS) 
corresponds to the time interval At of the correspond- 
ing Landau-Lifshitz-Gilbert equation in the high damp- 
ing limit |]l2|. Throughout the paper we set At =_Jq^s/i 
and adjust R, except of the simulations for Fig. [l] where 
we have to fix R and calculate At. 

In order to compute the dipole-dipole interaction effi- 
ciently we use FFT method s [jl7| f or the calculation of 
the long-range interactions [|l 4|, (l5|. This method uses 
the discrete convolution theorem for the calculation of 
the dipolar field 



H, 



E 



3ejj(e,y • Sj) Sj 



This dipolar field can be rewritten in the form 



(3) 



(4) 



where the Greek symbols r],8 denote the Cartesian com- 
ponents x, y, z and the Latin symbols i, j denote the lat- 
tice sites. (W v )ij are interaction matrices which only 
depend on the lattice. Since the lattice is translational 
invariant one can apply the discrete convolution theorem 
in order to calculate the Fourier transform of the dipolar 
field as 



(5) 



Thus one first computes the interaction matrices 
(W ve )ij and its Fourier transform {W v9 )k- This task has 
to be performed only once before the simulation starts 
since the interaction matrices depend only on the lattice 



structure and, hence, remain constant during the simula- 
tion. For each given spin configuration the dipolar field 
can then be calculated by first performing the Fourier 
transform of the Si, second calculating the product above 
following Eq. ||, and third transforming the fields Hfc back 
into real space resulting in the dipolar fields H^. Using 
FFT techniques the algorithm needs only of the order 
of NlogN calculations Jl7j instead of TV 2 calculations 
which one would need for the straightforward calculation 
of the double sum in Eq. [|. 

In order to apply the convolution theorem the system 
has to be periodically and the range of the interaction has 
to be identical to the system size. Since we are interested 
in finite systems with open boundary conditions we use 
the zero padding method |l^| . This method is based on a 
duplication of the system size where the original system 
is wrapped up with zero spins. 

The FFT methods is well established as containing 
no further approximation in the context of micromag- 
netism Jl8[ . There is however an additional problem 
concerning the implementation of this method to Monte 
Carlo algorithms since here the spins are not updated 
in parallel. Principally, in a Monte Carlo simulation the 
dipolar field has to be calculated after each single spin flip 
since the dipolar field at any site of the lattice depends on 
all other magnetic moments. Thus if one magnetic mo- 
ment moves, the value of the whole dipolar field changes. 
Hence, in a usual Monte Carlo algorithm one would store 
the dipolar fields in an array which is updated after ev- 
ery accepted spin flip |19|. Here, only the changes of the 
dipolar field due to the single spin update are calculated 
(N operations) resulting in a need of computation time 
which scales with N 2 per MCS. 
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FIG. 1. Reduced characteristic time T"//[i s vs. update in- 
terval At u of the dipolar fields for different trial step widths R. 
The solid line represents the "correct" value when the dipolar 
field is calculated after every spin update. The vertical lines 
represent the number of MCS which one spin needs at least 
to reverse. Dashed lines connecting the points are guides to 
the eye. L = 8, a = 1, k B T = 0.007J, /i s B 2 = 0.15J, and 
w = 0.032J. 
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Nevertheless, in the following we will show that alter- 
natively one can also recalculate the whole dipolar field 
at once after a certain number of MCS taking advantage 
from the FFT method. This can be a good approxima- 
tion as long as the changes of the system configuration 
during this update interval are small enough. 

In order to investigate the validity of this idea system- 
atically, we consider thermally activated magnetization 
reversal in a spin chain, in other words, in a cylinder 
with diameter D = 1. First, we calculate the "correct" 
value of the reduced characteristic time following from a 
Monte Carlo simulation where the dipolar field is calcu- 
lated after each accepted spin flip. This value is shown as 
solid line in Fig. [y. The data points represent the depen- 
dence of the reduced characteristic time on the length of 
the update interval At u after which the dipolar field is 
calculated. We used two different trial step widths R for 
comparison. 

As is demonstrated our method converges already for 
update intervals clearly longer than 1MCS, depending 
on the trial step width R. The dependence on R can 
be understood as follows. The vertical lines represent 
the minimal number of MCS which one spin needs to re- 
verse. This can be estimated from the mean step width 
of a magnetic moment jl2| assuming that each trial step 
is accepted following a strong external field. The mean 
step width within our Monte Carlo procedure is R/y/E 
and thus the minimal number of MCS for a spin rever- 
sal is Try/h/R. We conclude from this considerations that 
it is a good approximation to update the whole dipolar 
field after a certain interval of the Monte Carlo proce- 
dure which has to be smaller than the minimum number 
of MCS needed for a spin reversal. Note that it follows, 
that this method should never work for an Ising system. 
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FIG. 2. Computation time t vs. system size L of a spin 
chain. The data are from Monte Carlo simulations where the 
dipolar field is calculated either directly as explained in the 
text (o) or with FFT methods (A). 

Now we will turn to an investigation of the efficiency 
and capability characteristics of the method. A com- 
parison of the computation time needed when the dipo- 
lar field is cither calculated straightforward by updating 



the field after each accepted spin flip or with the FFT 
method after each MCS is shown in Fig. ||. Here, the 
computation time is the CPU time needed on an IBM 
RS6000/590 workstation for 100MCS. Using the FFT 
method the computation time is proportional to LlnL 
while for the usual method the computer time scales like 
L' 2 . For the chain of length 2 18 = 262144 the gain of effi- 
ciency of the FFT method is roughly a factor 5000 (note 
that the FFT algorithm is most efficient for the system 
sizes which are products of small prime numbers) which 
is a rather impressing result. 

As a test of the FFT method as well as the time 
quantification of the Monte Carlo algorithm we use 
for comparison also another numerical method namely 
Langevin dynamics simulations. Here we solve numer- 
ically the Landau-Lifshitz-Gilbert equation of motion 
with Langevin dynamics using the Heun method § . This 
equation has the form 



(1 + a 2 )n s dSj 
dt 
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= -S ( x H l (t) + a(S J x Hi(t))) (6) 



§£, the gyromag- 



with the internal field H;(t) = 
nctic ratio 7 = 1.76 x 10 n (Ts) -1 and the dimension- 
less damping constant a. The noise C >i represents ther- 
mal fluctuations, with (C;0)) = and (C?(*)C| (<0) = 
2SijS v eS(t — t')akBT^i s /j where i,j denote once again 
lattice sites and 77, 9 Cartesian components. The imple- 
mentation of the FFT method is straightforward for a 
Langevin dynamics simulation since here a parallel up- 
date of spin configurations is appropriate. 
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FIG. 3. Reduced characteristic time r"//ns vs. damping 
constant a. The data are from Monte Carlo and Langevin dy- 
namics simulations for a spin chain of length L — 64. The line 
resembles the (l+a 2 )/a behavior. jj, s B z = 0.15J, w — 0.032J 
and k B T = 0.025 J. 

In order to compare our different methods, in Fig. || 
the a dependence of the characteristic time for the re- 
versal of a spin chain is shown. Monte Carlo data are 
compared with those from Langevin dynamics simula- 
tions. Interestingly, for the whole range of a values the 
Monte Carlo and Langevin dynamics data coincide. This 
is in contrast to earlier tests jl2 ,|ll| where this agreement 
was achieved only in the high di 



amping limit. Seemingly, 
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there exist certain systems which show only a simple a 
dependence of the form r ~ (1 + a 2 ) /a. This a depen- 
dence usually describes the high damping limit which is 
rendered by the time quantified Monte Carlo algorithm. 
In our model this a dependence is obviously valid in the 
whole range of damping constants which we considered, 
leading to the remarkable agreement. In general, this is 
not necessarily the case since for other models different 
low and high damping limits appear (see e. g. |2(i|, f2l| ) , 
depending on the systems properties. In the following we 
set a = 1. 



100000 
10000 



1 1 1 






; 




MC O 


i i i 


LD • 

i i 



2 4 6 8 10 12 

AE/(k B T) 



IV. NUCLEATION AND CURLING 



As application of the methods described in the pre- 
vious section we start with a simple chain of classical 
magnetic moments of length L = 128 as extreme case 
of a cylindrical system with D = 1. Thermally acti- 
vated magnetization reversal in a spin chain was treated 
analytically by Braun |^2|, H) who proposed a soliton- 
antisoliton nucleation process as reversal mechanism. For 
a finite system with open boundary conditions the nu- 
cleation will originate at the sample ends leading to an 
energy barrier |24] 



AE = 2V2d(tanhr - hr) 



(7) 



with h — fi s B z /(2d) and r — arcosh-^/ \/h. Here, the 
dipole-dipole interaction is approximated by a uniaxial 
anisotropy of the form — d'^2, i {Sf) 2 with w w d/n [ pi] . 
This estimate follows from a comparison of the stray field 
energy of an elongated ellipsoid p8[ with the energy of 
a chain with uniaxial anisotropy. In |l3[| the occurrence 
of soliton-antisoliton nucleation was confirmed and Eq.|^ 
was verified as well as the asymptotic form of the corre- 
sponding characteristic time 



T exp 



AE 



(8) 



where tq is a prefactor depending on the system param- 
eters, the applied magnetic field, and the damping con- 
stant. This prefactor was also derived asymptotically for 
a periodic chain for the case of nucleation |^3| as well as 
for coherent rotation J25|. The latter is energetically fa- 
vorable for low enough system size. However, for a finite 
chain the prefactor is still unknown. 



FIG. 4. Reduced characteristic time t"//^i s vs. inverse tem- 
perature AE/ksT from Monte Carlo and Langevin dynamics 
simulations. The solid line is a fit to the low-temperature data 
in order to obtain the energy barrier numerically. L = 128, 
HsB z = 0.15J, and w = 0.032J 

Fig. [| shows the temperature dependence of the re- 
duced characteristic time for our spin chain where the 
dipole-dipole interaction is not approximated by a uni- 
axial anisotropy as discussed before jl3| but taken into 
account using FFT methods. Data from Monte Carlo 
simulations are shown as well as from Langevin dynam- 
ics. Both methods yield identical results. The slope of 
the curve in the low temperature limit which represents 
the energy barrier is 12% lower than the predicted one. 
This deviation is probably due to the local approximation 
of the dipole-dipole interaction underlying Eq. [?| 

We will now turn to the question whether such a nu- 
cleation mode will also be found in an extended system. 
Therefore we consider cylindrical systems of length 128 
with varying diameter D. Fig. ^| shows snapshots of two 
simulated systems at the characteristic time. For D = 8 
(left side) two nuclei originate at the sample ends lead- 
ing to two domain walls which pass the system. This 
is comparable to the nucleation process in a spin chain 
since all spins in a plane of the cylinder are more or less 
parallel (also shown in the figure) and can be describe 
as one effective magnetic moment leading to the above 
mentioned model of Braun (24) . 

The behavior changes for larger diameter above the so- 
called exchange length l ex . For a reversal mode with in- 
homogeneous magnetization within the planes a domain 
wall has to be created. The loss of energy due to the 
existence of a 180° domain wall on the length scale I in 
a spin chain is 



AE = J^(l - cos(6>; - 6jj) 



Jir 2 
21 



(9) 



under the assumption that the change of the angle 6 be- 
tween the next nearest magnetic moments is constant (in 
a continuum description one can also prove that this re- 
sembles the minimum energy configuration by minimiz- 
ing the Euler Lagrange equations). 
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FIG. 5. Snapshots of simulated spin systems at the characteristic time for nucleation (left) and curling (right). The z com- 
ponent of the magnetic moments is also color coded (green corresponds to spin down, blue to spin up). The clipping plane on 
the left side shows a layer in the upper domain wall that on the right side shows a layer with a magnetization vortex in the 
central curling region. y, s B z = 0.15J, w = 0.032J, L = 128. D = 8 (left) and 16 (right). 
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The dipolar field energy of a chain of length I is at 
most 3iuZ£(3) where we used Riemann's Zeta function 



oo 1 



(10) 



(see also |l!| for a corresponding calculation in two di- 
mensions) . A comparison of Eq. ^| and Eq. [fT^ yields the 
exchange length 



J 



6u<(3) 



(11) 



Usually |18,24 the exchange length is calculated in the 
continuum limit where 3£(3) ~ 3.6 is replaced by 7r. How- 
ever, we prefer the slightly deviating expression above 
since its derivation is closer to the spin model which we 
discuss here. The exchange length for our material pa- 
rameters is l ex « 7, so that for systems with smaller 
diameter it cannot be energetically favorable for the sys- 
tem to build up inhomogeneous magnetization distribu- 
tions within the planes while for wider systems another 
reversal process can occur namely curling |26| , |27f| . 

A snapshot during the reversal process in a wider par- 
ticle (D = 16) is also shown in Fig. ||. The middle part 
of the cylinder is in the curling mode consisting of a 
magnetization vortex with a central axis which is still 
magnetized up. The reversal mode shown in this figure 
is still mixed up with an additional nucleation process 
which started first at the sample ends. Note, that in the 
clipping plane on the right hand side one can find the ex- 
change length as typical length scale of the domain wall 
in the vortex. 



V. CONCLUSIONS 

We considered thermal activation in classical spin sys- 
tems using Monte Carlo methods with a time quanti- 
fied algorithm as well as Langevin dynamics simulations. 
We combined both techniques with FFT methods for the 
calculation of the dipolar field. Whereas this method 
is established in the context of micromagnetic equations 
of motion it is less obviously applicable to Monte Carlo 
methods with single spin flip dynamics. We show that it 
can be a good approximation for vector spins to update 
dipolar fields during a Monte Carlo procedure only after 
certain intervals the length of which depend on the trial 
step width of the algorithm. Since Monte Carlo methods 
need less computational effort as compared to Langevin 
dynamics simulations its combination with FFT meth- 
ods strongly enhances the capabilities for the numerical 
investigation of thermal activation. 

As application we consider a model for nanowires. 
We compare numerical results for the characteristic time 
of the magnetization reversal in a spin chain with the 



theoretical formula for the energy barrier of soliton- 
antisoliton nucleation. We found small deviations proba- 
bly due to the fact that Braun in his model approximated 
the dipole-dipole interaction by a uniaxial anisotropy 
pij . Varying the diameter in a cylindrical system we ob- 
serve a crossover from nucleation to a curling-like mode. 
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